function [derivs] = del_log_det_B(hyps, func, n, n_class, X, y, approxF)
% a function evaluating d log.det B d theta_i
% by Mark Norrish, 2011

% hyps: as a vector; [ h_class1 ; ... ; h_class_c ]

dim = length(hyps) / n_class;
Hyps = reshape(hyps, dim, n_class);

bigK = zeros(n*n_class); K = zeros(n,n,n_class); sigma_noise = 1e-7;
for c = 1:n_class
  K(:,:,c) = func(Hyps(:,c), X, X) + sigma_noise*eye(n);
  bigK(1+(c-1)*n:c*n,1+(c-1)*n:c*n) = K(:,:,c);
end

if nargin <= 6
  f = alg_3_3(n, n_class, K, y);
else
  f = alg_3_3(n, n_class, K, y, approxF);
end
%f= approxF;
F = reshape(f, n, n_class); % it's SO much easier to think with
expsum = repmat(sum(exp(F)')',n_class,1);
pi = exp(f)./expsum;
bigPi = zeros(n_class*n,n);
for i = 1:n_class
  bigPi((i-1)*n+1:i*n,:) = diag(pi((i-1)*n+1:i*n));
end
W = diag(pi) - bigPi * bigPi';
Pi = reshape(pi,n,n_class);
derivs = zeros(size(hyps));

expsum = expsum(1:n);
dpyfdf = reshape(y, n, n_class) - Pi; % d log p(y|f) / d f
IKW = (1 + sigma_noise) * eye(n*n_class) + bigK * W; % this is aka B
L = chol(bigK) + chol(W + inv(bigK));% IKW^0.5; % now L'*L = IKW
del_logq = zeros(1,n); % dlog(q)/d(fhat_i^c)
for c = 1:n_class
  crn = 1+(c-1)*n:c*n; % class-range: i.e. the class-appropriate block
  expf = exp(F(:,c));
  for i = 1:n
    equalities = repmat(c == 1:n_class,n_class,1);
%    del_W = zeros(n_class);
%    for c1 = 1:n_class
%	for c2 = 1:n_class
%	  del_W(c1,c2) = -Pi(i,c1)*Pi(i,c)*(c1 == c2) + 2 * Pi(i,c)*Pi(i,c1)*Pi(i,c2) - ((c==c1) + (c==c2)) * Pi(i,c1)*Pi(i,c2);
%	end
%    end
    del_W = 2 * Pi(i,c) * Pi(i,:)' * Pi(i,:) - (equalities + equalities') .* (Pi(i,:)'*Pi(i,:)) - diag(Pi(i,:)*Pi(i,c)); % dW / df
%    fprintf('Error%f\n', sum(sum(abs(my_del_W - del_W))));
    del_W(c,c) = del_W(c,c) + Pi(i,c);
    z = zeros(n); z(i,i) = 1;
    del_W = kron(del_W, z);
    del_logq(i) = trace((L \ (L' \ bigK)) * del_W);
  end
  for hyp_idx = 1:dim
    h = (c-1)*dim+hyp_idx;
    del_K = zeros(n*n_class);
    del_K(crn,crn) = func(Hyps(:, c), X, X, hyp_idx);
    del_f = (L \ (L' \ del_K))(crn,crn) * dpyfdf(:, c); % d fhat / d theta_i; by eq 5.24
    % eq 5.22: gimped down 
    derivs(h) = -0.5 * trace((L \ (L' \ del_K)) * W) - 0.5 * del_logq * del_f; % explicit, implicit
  end
end

derivs = -derivs; % because it's negloglik not loglik